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Widely-used forms of the K-e turbulence model are shown to yield arbitrary steady- 
state converged solutions that are highly dependent on numerical considerations such 
as initial conditions and solution procedure. These solutions contain pseudo-laminar 
regions of varying size. By applying a nullcline analysis to the equation set, it is pos- 
sible to clearly demonstrate the reasons for the anomalous behavior. In summary, the 
degenerate solution acts as a stable fixed point under certain conditions, causing the 
numerical method to converge there. The analysis also suggests a methodology for 
preventing the anomalous behavior in steady-state computations. 

I. Introduction 

The prediction of turbulent flow fields for engineering purposes continues to be dominated 
by employing the Reynolds-averaged Navier Stokes (RANS) approach. Although there are sev- 
eral classes of closure methodologies available, the class of two-equation linear eddy viscosity 
models may be the most popular. Included in this class is the K-e model, which has many 
variants. Such models are most effectively utilized when the focus is on mean field dynamics 
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rather than detailed behavior of the statistical moments. The two-equation formulation yields an 
eddy viscosity that directly influences the mean flow behavior. The K equation can be directly 
derived from the full Reynolds stress transport equation (by taking the trace). Closure of this 
K equation then requires a constitutive equation for the Reynolds stress tensor as well as the 
transport and pressure-diffusion terms. For the Reynolds stress tensor, a linear relation with the 
mean strain rate tensor is assumed with the proportionality coefficent being the turbulent eddy 
viscosity. For the transport and pressure-diffusion terms gradient transport models are generally 
assumed. The e equation is based on transport processes in the dissipative range, but it is often 
viewed as being an empirical model . 1 - 2 

Naturally, implicit in the use of all RANS models is the assumption that for a given set of 
initial and boundary conditions, a unique solution will be obtained. However, it will be shown 
here that for certain formulations of the K-e model, portions of the flow field can converge to 
a degenerate solution that is “pseudo-laminar” in nature. The terminology “pseudo-laminar” 
alludes to the fact that the model does not predict the correct laminar limit in terms of the tur- 
bulent to mean flow time scale ( SK/e ), but the resulting eddy viscosity is still near zero so the 
mean flow behaves as a laminar flow. Unfortunatley, this disturbing behavior is further exas- 
cerbated by the fact that the location and spatial extent of these pseudo-laminar regions in the 
converged solutions can depend on initial conditions and method of solution! This paper both 
demonstrates the problem and performs an analysis that explains why the system of equations 
behaves in this manner. It should be noted that this anomalous behavior occurs only in develop- 
ing flows, i.e., boundary layers. Flow configurations that utilize streamwise periodic boundary 
conditions (such as fully developed channel flow) do not in general encounter this problem. 

Such issues have not been addressed previously, and the purpose of this analysis is to criti- 
cally examine the characteristics of the K-e model from a dynamical systems standpoint. Mo- 
hammadi and Pironneau 3 reported extensive mathematical analysis on the K-e model, but did 
not consider the anomalies described here. Dynamical system analyses have previously been 
utilized 4-6 in order to identify the fixed points of two-equation and second-moment closures 
in homogeneous shear flow, and to calibrate closure models for equilibrium flows. The issues 
addressed in the present study, however, require that the evolution toward the equilibrium state 
is understood, particularly the solution trajectories from a given initial condition, and not the 
fixed points per se. Although the theoretical approach is inherently based on the simplifying 
assumption of homogeneous turbulence, an attempt is also made to account for inhomogeneous 
effects so that the theoretical results can be more easily related to the full numerical solution of 
the model in realistic wall-bounded flows. 
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As already alluded to, the dynamical systems analysis undertaken here is not only con- 
cerned with fixed point identification but also about dependent variable phase-plane trajectories 
and dependence on initial conditions. The analysis is shown to isolate deficiencies in some 
formulations of the K-e model that lead to initial condition and solution-method dependent 
converged solutions. 


II. Illustration of the Problem 


To illustrate the anomalies arising in numerical solutions of the K-e model, two cases are 
considered: (1) a flat plate boundary layer flow, and (2) flow over an airfoil. The computer 
code CFL3D 7 was employed. 3 As will be shown, the capacity to reach arbitrary steady-state 
solutions is a property of the K-e equations themselves, so any numerical method will encounter 
the problem. 

A basic form of the K-e model can be written as 
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where D/Dt = d/dt + Ujd/dxj , C e i = 1.44, C e2 = 1-83, cr K = 1.0, a e = k 2 /[y^(C s2 ~ 
Cei)], k = 0.41, and = 0.09. The set of values used here for the closure constants are 
representative of values typically chosen. However, as will be seen, the anomalous behavior 
exists regardless of their exact numerical values. The rate of turbulent energy production V is 
defined as —uiuj dUi/dxj, and the components of the Reynolds stress tensor is modeled using 
the Boussinesq assumption, 


UiUj = —I\5jj — v t 
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with the eddy viscosity given by 

r2 

Vt = — . (3) 

e 

In this type of formulation, the function f 2 is introduced into Eq. (lb) so that as a solid 
boundary is approached the destruction-of-dissipation term does not become ill-behaved. 9 It is 

a Another code, ISAAC, 8 has also been used and yields similar results. 
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also used to calibrate the model in the log-layer region. For example, two of the most widely- 
used forms for f 2 are: 

.I 2 = I - ci exp (— c 2 Re \ ) , (4) 

where Re T = K 2 /(ve), and 

f 2 = 1 - c 3 exp ( -c 4 Re K ) , (5) 

where Rex = \Z~Kdjv and d is the distance from the wall. The constants ci and c 3 are always 
positive and less than or equal to 1.0, and c 2 and c 4 are positive constants. Note that 1 — c\ < 
f 2 < 1 for Eq. (4) and 1 — c 3 < / 2 < 1 for Eq. (5). For the present demonstrations, only Eq. (5) 
was used, with the commonly used values c 3 = 1, c 4 = 2/25, but Eq. (4) exhibits similar 
behavior. 

In order to rule out the possibility that the observed arbitrary behavior occurred only for 
a particular set of boundary conditions, various freestream boundary conditions for K and e 
were employed. The anomalies occurred regardless of these variations, as exemplified below. 
Because the solutions were highly dependent on numerical parameters, the grid size itself could 
also have an influence on the final solution. Reasonably fine grids were used for both cases, 
but a formal grid independence study was not conducted because it has no meaning when the 
equations themselves (and not the numerics) can yield arbitrary results. 

A. Flat Plate Boundary Layer 

In the flat plate boundary layer flow considered first, a freestream turbulent intensity of 0.2% 
(= y/2K 00 j‘i) was imposed everywhere in the computational domain using a grid of 193 x 65 at 
a Reynolds number (based on plate length) of Re = 6 x 10 6 . Thus, the initial condition K 0 on K 
was everywhere the same as the boundary condition (7% = l\' y = l\\ y Ju 2 y = 6.0 x 10 -6 ). The 
dissipation rate boundary condition at inflow, which determines the freestream eddy viscosity 
and turbulence freestream decay rate, was set at l/u^ = 8. 1 /'«( c , or e' x: = 8. 1 /\/ c . This 

yielded a freestream eddy viscosity for this case of //,/ //,- x: = 0.40 at the boundaries. A sample 
solution showing the final u/U ^ -velocity contours is plotted in Fig. 1, for the case with initial 
condition e[, = 8.1 7%. The solution exhibited a pseudo-laminar solution upstream of x /l « 0.1, 
then a turbulent solution downstream. Next, different cases were run with four different initial 
conditions % in the field varying from 0.00817% to 81/%. 

As Fig. 2 shows, each converged solution yielded a different apparent “transition” location 
that was located further downstream with increasing initial dissipation rate values. At suffi- 
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ciently high levels (£q = 81 1\' 0 ), the flow remained laminar throughout. Initially, this might 
seem like consistent behavior: an increased dissipation should reduce the level of turbulent 
kinetic energy K and therefore should shift the position of “transition” further downstream. 
However, this rationale is flawed for two reasons; one more serious than the other from a CFD 
practitioner’s point of view. 

First, in the presence of mean shear (e.g. S = dU/dy), a laminar solution is characterised 
by very large values of the time-scale ratio, i.e., SK/e >■ 1. Physically this implies that the 
turbulent time scale is much greater than the mean flow time scale and as such turbulence does 
not persist. But inspection of the “laminar” regions in the current solution shows that SK/e 
is in fact < 1. In other words, the I\-e model is not really predicting laminar flow at all, but 
rather a pseudo-laminar behavior that will be shown in section III to be a stable fixed point of 
the equations. 

Second, and even more serious is that the final converged solution should not depend on the 
initial condition in a steady-state computation at all! The converged solution should depend 
on the boundary conditions only, and in this case the boundary conditions were the same in 
all computations. It should be stressed that the solutions presented in Fig. 2 were very well 
converged solutions. The /.2-norm of the density residual dropped by more than 6 orders of 
magnitude, to approximately 1 x 10 — 14 , as shown in Fig. 3. The solutions after 25,000 multigrid 
cycles showed no perceptible differences from those solutions obtained after 2500 cycles. 

B. RAE 2822 Airfoil 

As a second example of anomalous behavior, the flow over an RAE 2822 airfoil at freestream 
Mach number Moo = 0.75, angle-of-attacka = 2.72, and Re = 6.2 x 10 6 was computed. A plot 
showing the airfoil shape and resulting pressure contours for these conditions is given in Fig. 4. 
There is a strong shock wave present on the airfoil upper surface near 65% chord, whereas the 
flow on the lower surface remains subsonic. In this example, the initial and boundary conditions 
were kept fixed but the numerical solution method was changed. The initial conditions and 
farfield boundary conditions in this case were set to l\' x — 2.5 x 10 -8 and e F = 3.75 x 10 -8 . 
This corresponded to a very low freestream turbulent intensity of 0.013% and a freestream eddy 
viscosity of p t /poo = 0.009. These are typical values used in CFL3D. 7 

Figure 5 shows the skin-friction distribution on the lower surface of the airfoil obtained using 
two different numerical solution strategies for obtaining converged solutions. The converged 
results were completely different, with each suggesting a “transition” location in a different 
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place. The first case was run using multigrid and 3 levels of mesh sequencing, with 2500 
iterations on the coarse grid, followed by 2500 iterations on the medium grid, and finally 3000 
iterations on the finest grid (257 x 97). The second case was run with multigrid and 2 levels of 
mesh sequencing, with 5000 iterations on the medium grid followed by 3000 iterations on the 
finest grid. Although not shown, both cases converged very well, with the /. 2 -norm of density 
residual reduced more than 3 orders of magnitude. 

Both these examples show that caution needs to be exercised when using the K-e model. 
A numerically converged solution does not necessarily constitute the intended solution to the 
set of governing equations; it may depend on numerical parameters such as initial conditions 
and solution procedure. It is also important to mention here that many CFD practitioners have 
noticed that the K-e equations often fail to go fully turbulent, although the cause has never 
been identified before. In fact, it is customary to build in ad-hoc fixes to attempt to ensure that 
turbulence always develops. Some of these fixes include: (1) restarting K-e solutions from 
another turbulent solution, (2) setting initial conditions to have turbulent-like levels rather than 
freestream levels, and (3) imposing a temporary source term in the boundary layer. All of these 
fixes, in general, are workable ways to avoid the problem; but they do not shed any light on 
the reasons behind the problem and were not developed based on any firm rational foundation. 
As a consequence, their generality cannot be assured. In the following section an analysis is 
conducted that makes the reasons clear. 


III. Dynamical Systems Analysis 


A dynamical systems analysis can be used to determine the temporal dynamics associated 
with the numerical solution of systems of equations. 10 A so-called nullcline analysis will also 
be used to identify some parametric restrictions on the K-e equations Eq. (1), in order to avoid 
arbitrary pseudo-laminar converged solutions. 

It is possible to gain critical insight into the solution behavior for inhomogeneous turbulent 
flows through an analysis of the homogeneous form of Eq. (1). In its homogeneous limit, Eq. (1) 
can be written as a nonlinear, autonomous equation system in the form 


dK* 

dt* 

de 

dt* 


c, 


I\ 


*2 


— e 


C el C,K* - f 2 c s2 — , 


(6a) 

(6b) 
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where K* = SK ( S = du/dy is fixed in the analysis), and t* — St. There can be either one 
or two critical points in this system: (;) the null vector with elements I\* = e = 0, and (ii), if 
it exists, is the intersection of the set of points with dK* jdt* = 0 that lie on the /\'*-nullclinc 
described by 

e = ±y/U^K*, (7a) 

and the set of points with dejdt* = 0 that lie on the e-nullcline described by 



Realizability considerations dictate that only the positive roots need to be considered. At the 
intersection point of Eq. (7a) and Eq. (7b), f 2 = C £l / C e2 . Thus, this second critical point 
exists only if f 2 can achieve a value C el / C s2 (which for the current set of closure coefficients is 
0.7869). From Eq. (4) or Eq. (5), this means that the second critical point (ii) can exist only if 
ci > 1 — C e i/C e2 or c 3 > 1 — C el /C e2 for these particular choices of f 2 . 

If stable, the critical points represent the possible steady- state solutions to the system, 
Eq. (6). Note that neither of these critical points is the so-called “turbulent” solution. In the anal- 
ysis of homogenous turbulent flow, the “turbulent” solution grows without bound ( K * — > oo). 
The reason the analysis yields an unbounded growth in this case is that the mean flow field is 
fixed and unaffected by the turbulence, and this provides an infinite source of energy for the 
turbulence. In practical computations, this behavior is not seen because there is a two-way 
coupling between the turbulent and mean flow fields. The coupling allows for a steady-state 
to be reached. Diffusion, to be introduced later in the analysis, also plays an important role in 
practical computations. 

The stability properties of the two critical points can be examined by linearizing about each 
critical point. The coefficient matrix of this linear system is the Jacobian matrix 


2C„ (iR 

Cfj.C s i + f 2 C e2 — C e 2 


-c,(/F) 2 - 1 

—2f 2 C s2 ( T^r ) - t'e’l (^) 


To determine the nature of the critical points (e.g., if they are stable or unstable), the eigenvalues 
of this matrix are found at the critical points. Table 1 lists the possible types. A center indicates 
that trajectories orbit around the critical point. A stable critical point means that, when solving 
the equations, the particular point can be reached; an unstable critical point or a saddle point 
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means that a numerical scheme will not converge to it (a saddle point is not stable in prac- 
tice because orbits approach the critical point along one eigenvector, but then recede along the 
eigenvector associated with the unstable solution). 

Table 1. Type of critical point as determined by eigenvalues of J. 


type of eigenvalues 

critical point type 

complex with zero real part 

center 

complex with negative real part 

stable focus 

complex with positive real part 

unstable focus 

real and both negative 

stable node 

real and both positive 

unstable node 

real with one positive, one negative 

saddle point 


It can be shown that, when it exists, the second critical point (z'z) of the system Eq. (6) is 
always a saddle point. When Eq. (4) is used for f 2 , the eigenvectors at the saddle point are: 


X x 


2 

X 2 


V / C^( 1 + c el ) + C e2 -s/U^D ln( A) ± y/A 


where ,4 = C,( 1 - C sl -) 2 + C,C € 2 D 2 W(B) + 8 G p C#D{C el - 3) In (B), B = (C s2 - 
C s i)/(cxC e 2 ), and I) = 1 — CA / C e2 • When Eq. (5) is used for / 2 , the eigenvectors at the 
saddle point are: 


X x 


2 

x 2 


V^(l + Csi) ± ^C,{l-C £i y-C,C s2 D\n(E) 


( 10 ) 


where E = (C e2 — C e i)/(c 3 C e2 ). The effect of the saddle point on the solution trajectories 
will be shown below when phase plots are drawn. But the more interesting critical point in this 
analysis is the first (z) degenerate one (K* — e — 0), which corresponds to a pseudo-laminar 
solution. When Eq. (4) or Eq. (5) is used for f 2 , and if the following condition is true, 

K* 

< 1 as A* = e = 0 is approached, (11) 


then this degenerate point is either a center or a stable point. This stability of the degenerate 
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critical point turns out to be the cause of the apparently arbitrary solution behavior demonstrated 
in section II. 

In section C, the conditions for which Eq. (11) occur will be explored in greater detail. It will 
also be shown that stability of the degenerate point requires the second critical (saddle) point 
to exist. For now, however, the assertion is made that in practice Eq. (11) is quite often true 
for the system of equations given by Eq. (1) or Eq. ( 6 ). As a result, Rey <§C 1 and Rex <§C 1 
near the critical point as well. Therefore, for ^ 1 Eq. (4) is approximately f 2 = 1 — c\, 
and for c-, — 1 it is approximately f 2 = c 2 K* 4 /( S 4 v 2 e 2 ). Similarly, for c 3 ^ 1 Eq. (5) is 
approximately f 2 = 1 — c 3 , and for c 3 = 1 it is approximately f 2 = c 4 l\ x] l 2 dj (uS^ 2 ). Using 
these expressions, along with expressions for the derivatives of Eq. (4) and Eq. (5) with respect 
to K* and e, the eigenvalues of the matrix in Eq. ( 8 ) near the degenerate critical point can be 
determined. These are given in Table 2 for the two f 2 expressions and various combinations of 
their coefficients. In all cases but one, the degenerate critical point is stable. (The “center” type 
of critical point is also not desirable in this context, because solutions can become locked in an 
orbit and fail to converge. However, no known models actually use c\ — 1 in Eq. (4), so from 
now on the analysis focuses only on the more common cases where the critical point is stable.) 

Table 2. Eigenvalues near the K* =6 = 0 critical point, with K*/e 1. The computations in section II 
employed a model that corresponds to the third alternative. 


f 2 function (c 2 , c 4 > 0 ) 

type of eigenvalues 

critical point type 

1 — Ciexp(— c 2 Rej ), 

Ci = 1 

complex, near- zero real part 

center 

1 — C!exp(— c 2 Re\), 

0 < ci < 1 

real, both negative for ci < (C e2 — 1 ) / C e2 
complex, negative real part for ci > (C s2 — 1 )/C s2 

stable node 
stable focus 

1 - e 3 (’x p( c 4 R< k), 
c 3 = 1 

complex, negative real part 

stable focus 

1 — e 3 ex|)( c 4 /iV/ v -), 
0 < c 3 < 1 

real, both negative for c 3 < (C s2 — 1)!C e2 
complex, negative real part for c 3 > (C s2 — 1 ) / 6 V 2 

stable node 
stable focus 


Equations (7a) and (7b) can be sketched for a case when the two nullclines intersect. In- 
spection of Eq. ( 6 a) and Eq. ( 6 b) then shows the following to apply: 


dK 1 

dt* 


= 0, on Eq. (7c 


dK 1 

dt* 


< 0, left of Eq. (7a) 


dK 
dt * 


> 0, right of Eq. (7a) (12a) 
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de 

dt* 


0, on Eq. (7b) 


- - — < 0, above Eq. (7b) 


> 0, below Eq. (7b). (12b) 


Combining this information with the knowledge of the behavior at the critical points, a clear 
mapping of the phase-plane trajectories can be drawn, as shown in Fig. 6. The eigenvectors 
(Xi, X 2 ) at the saddle point are shown along with dividing curves (dashed lines) that separate 
regions of different trajectory behavior. Clearly, any initial condition above dividing curve 1 
will eventually end up at the degenerate critical point. 

An actual phase-plane portrait can be constructed by computing the right hand sides of 
Eq. (6) for a large number of K* and s values. These computed values then correspond to 
dK* /dt* and de/dt* at each particular point in nondimensional K*-e phase space, and the 
trajectory (how I\* and e change with time or with iteration) can be computed as well. 

An example is shown in Fig. 7, using Eq. (5) for f 2 with c 3 = 1, c 4 = 2/25, S' = 1569, 
d' = 2.945 x 10 -4 (typical results using Eq. (4) yield a similar behavior). The nullclines are 
shown along with the phase space trajectories. It is clear that this figure matches the sketch 
in Fig. 6 in character: K* = e = 0 is a stable attractor, and the other critical point (near 
I\* = 0.19, e = 0.057 in this particular case) is a saddle point. The true turbulent solution is 
obtained when K* grows (exponentially); in other words, the expected turbulent solution only 
occurs when the K* and e values follow the trajectories in the lower right or far upper right 
parts of the plot. This figure demonstrates that there are many regions in the map for which the 
solution converges toward the degenerate critical point K* — e — 0. 

It is also interesting to look at the phase space trajectories of the case with f 2 — 1, shown in 
Fig. 8. Here, the two nullclines do not cross, so there is no saddle point (i.e., no second critical 
point). The trajectories now behave like the ones to the right of the saddle point in Fig. 7; 
regardless of the initial condition, the solution always goes to the turbulent solution and never 
to the degenerate one. 

A. Accounting for Diffusion 

Inhomogeneous effects that necessarily occur in any practical calculation can be represented 
by the viscous diffusion term and the gradient diffusion models for the turbulent transport, 
represented by vt/&K and v t /a e . Since the consideration here is for high Reynolds number 
flows, viscous diffusion effects can be neglected in the turbulent region, and it is only necessary 
to adequately represent the effects of the turbulent tranport in both the kinetic energy and energy 
dissipation rate equations. It suffices for the purposes here to simply assume that the transport 
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effects act over a distance l given by /\' 3,/2 /£. Thus, qualitative estimates for the transport terms 
can be written as 


d ( v t dK 


dx j \ a k dxj 
d ( v t de 


dxj 




dx. 




O(e) 

• = 

£ 2 

(13a) 

k) 

,l£ ~K ’ 

(13b) 


with C* K and C* e as unknown coefficients. Using these estimates for the transport terms al- 
lows them to be grouped with the destruction terms. Thus, an equation set that accounts for 
inhomogeneity can be written as 


dK * 
dt* 
de 
dt* 


K* 2 

C, (1 -c; K )e 

£ 

C'lW - (f 2 c e2 - c;^ . 


(14a) 

(14b) 


It is clear to see that the effect of the additional diffusion-type terms is merely to tilt the 
nullclines up, making them steeper. As will be shown in section B, this effect is needed 
to achieve agreement between theory and computation. The relative shapes and positions 
of the phase space trajectories are a function of the relative magnitudes of the ( a priori un- 
known) C* K and C * s . In this case, the nullclines intersect to form a saddle point when f 2 = 
[c.i(i-c; K )+c;,]/c.2. 


B. Comparison of Computations with Analysis 

Through the foregoing analysis, it is clear that the equations themselves (when they contain 
an f 2 function in the e-equation destruction term) can cause degenerate solutions to occur. But 
how well does the theory compare with actual RANS computations? Although the actual RANS 
computations are much more complicated than the analytical model because the mean flow S 
and the actual diffusion terms vary in time and space, we can look at these varying values at 
specific points from the earlier flat plate computation, for example, and choose representative 
levels over the latter part of the temporal development. These representative levels can then be 
inserted into Eq. (14) when computing the theoretical phase-space trajectories for comparison. 

Figures 9 and 10 show the results using this procedure at two different points in the flow 
field. The first figure shows computed results at a point in the boundary layer that eventually 
converged to a pseudo-laminar degenerate result, along with the theoretical trajectories. Here, 
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diffusion effects were negligible, and there is excellent agreement between theory and compu- 
tation. The second figure is for a point further downstream that became fully turbulent. Here, 
the effect of diffusion was important and needed to be taken into account in the theory; values 
used for C* K and G* e were 0.28 and 0.59, respectively, based on representative levels seen in 
the flat plate computation. There is again excellent agreement between theory and computation. 


C. Avoiding Arbitrary Steady-State Solutions 

Earlier, Eq. (11) was given as a condition for which the degenerate critical point was either a 
center or a stable point. The interrelationship between this condition and the value of f 2 is now 
explored. Consider the equation for I\* /e: 


d(K*/e ) 
dt* 


--((?,! - 1)-(1 -C. 2 / 2 ), 

£ 


(15) 


where Vje = C^K* 2 [e 2 . From this equation it is immediately apparent that for d( K */e)/dt* < 
0 (which must be true for K* /e to be driven toward zero), the following must hold: 


, VC S i-ll 

J2 < — — ^ r ■ 

£ W2 V'e2 


If this inequality is not satisfied, then K* je remains finite, and the eigenvalue analysis of the 
Jacobian matrix Eq. (8) shows that the K* — e = 0 critical point is unstable, and no longer an 
attractor for the degenerate solution. This suggests a method to avoid the problem of arbitrary 
solutions for steady-state computations: compute / 2 as usual with Eq. (4) or Eq. (5), but then 
limit it via 

„ . r, ( t vc el ~ i , i m 

J-2 = mm 1, max / 2 , — 77- , ( 17 ) 

V £ W2 C s2 y_ 

during the early transient stages of a steady-state computation. This limiter does not allow the 
value of / 2 to go below the critical level defined by Eq. (16). Once turbulence has been estab- 
lished, then the limiter can be removed. This analysis does not lead to a method for avoiding 
arbitrary solutions in time-dependent computations. Such computations are prone to anomalous 
behavior as well. 

Another important observation is that the right hand side in Eq. (16) is less than Cel / C S 2 
for V/e < 1. Because the second critical (saddle) point exists only if / 2 can achieve a value 
C £ i/C £ 2 , this implies that a second critical point is necessary in order for the degenerate point 
( K * = £ = 0) to be stable, when using the current equation set. 
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D. Extension to More General Form 


To summarize, the analysis suggests that the presence of two critical points for a commonly- 
used form of the K-e equations has the potential to produce anomalous pseudo-laminar behav- 
ior. The degenerate point associated with the pseudo-laminar solution (A * /e -C 1 ; /\ * = £ = 0) 
is stable only if there exists a second critical point given by f 2 = C s i/C e2 . 

One can write the homogeneous limit of the K-e equations in more general form: 


dK* 
dt* 
da 
dt * 


K* 2 

c: e 


■ & 
£ 

lx* 




K 


*2 



(18a) 

(18b) 


where the new variables C * , C* x , and C* 2 can now each include functions of the solution (for 
example, low-Reynolds-number models 9 often use a function [,, that is a part of C* ) . In the 
more general form Eq. (18), the second critical point is now defined by 


C*e2 = C* s V 


(19) 


The analysis strongly suggests that any I\-e model for which C* 2 A C* x somewhere in the flow 
has the potential to yield an arbitrary pseudo-laminar solution. 


IV. Concluding Remarks 

A peculiar problem inherent in a widely-used form of the K-e turbulence model has been 
demonstrated and analyzed. This problem has a potentially large impact on practical CFD com- 
putations. Use of an f 2 function multiplying the destruction term of the dissipation rate equation 
was shown to cause portions of the flow field to converge to a degenerate pseudo-laminar con- 
dition. Most disturbingly, this condition is highly dependent on numerical parameters such as 
initial conditions and solution procedure. In other words, RANS solutions using this particular 
form of the K-e model can easily yield arbitrary fully-converged solutions with pseudo-laminar 
regions of varying size. Time-dependent computations are also susceptible to the anomalous be- 
havior. 

A nullcline analysis was used to analyze the homogenous form of the equations, followed by 
a form that approximately accounts for the effect of diffusion. The analysis clearly demonstrated 
the reasons for the anomalous behavior of this turbulence model: the degenerate solution was 
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a stable fixed point under certain conditions, causing the numerical method to converge there. 
The analysis also led to a methodology for preventing the anomalous behavior in steady-state 
solutions using the current equation set. 

The results presented here also suggested that any K-e model for which the coefficient 
multiplying the destruction term in the e equation can be less than or equal to the coefficient 
multiplying the production term has the potential to produce arbitrary pseudo-laminar solutions. 
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Figure 1. Example flat plate solution showing u/lJ„_ contours for initial condition s' 0 = 8.1 K' 0 ; M,-, = 0.2, 
Re = 6.0 x 10 6 (I = length of plate). 
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Figure 2. Streamwise variation of skin-friction coefficient on front half of flat plate as a function of initial 
conditions (K' 0 — 6.0 x 10 -6 ). 
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log(L 2 -norm of residual) 



Figure 3. Convergence history for flat plate computations. 
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Figure 5. Streamwise variation of skin-friction coefficient on RAE 2822 airfoil lower surface for two differ- 
ent solution procedures (both procedures converged). 


19 of 24 





Figure 6. Sketch of nullclines and critical points from Eq. (7a) and Eq. (7b) and resulting K*-e phase 
diagram (K* = SK); dividing curves (dashed lines) divide trajectories near critical points. 
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Figure 7. Example phase-plane portrait of Eq. (6) showing nullclines (two oblique lines going from lower 
left to upper right) and trajectories (lines with arrows) with f 2 given by Eq. (5). 
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Figure 8. Example phase-plane portrait of Eq. (6) showing nullclines (two oblique lines going from lower 
left to upper right) and trajectories (lines with arrows) with f 2 — 1. 
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Figure 9. Comparison of theory with computed result at a point in pseudo-laminar region of the flat plate 
computation. 
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Figure 10. Comparison of theory with computed result at a point in turbulent region of the flat plate com- 
putation. 
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